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The universal finite-size scaling function of the critical Casimir force for the three dimensional XY 
universality class with Dirichlet boundary conditions is determined using Monte Carlo simulations. 
The results are in excellent agreement with recent experiments on 4 He Films at the superfluid 
transition and with available theoretical predictions. 
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Casimir forces are always present in nature when a 
medium with long-range fluctuations is confined to re- 
stricted geometries. The quantum mechanical Casimir 
effect was proposed theoretically 60 years ago by 
H. B. G. Casimir [l| and describes an attractive force be- 
tween two conducting plates in vacuum, induced by the 
vacuum fluctuations of the electromagnetic field. Fur- 
thermore, Goldstone modes 0] and surface fluctuations 
[|| can give rise to Casimir forces. Near continuous phase 
transitions, long-range fluctuations of the order parame- 
ter lead to the analogous thermodynamic Casimir effect 
4] , which can change the thickness of critical liquid films 
5(. In a series of papers, Garcia and Chan [f| Q and 
Ganshin et al. Q were able to measure the thinning of 
liquid 4 He films close to the A-point due to the thermody- 
namic Casimir effect. They found a characteristic deep 
minimum (dip) in the film thickness just below the super- 
fluid transition temperature T\. Using finite-size scaling 
methods, they accurately determined the scaling function 
•d{x) of the Casimir force, which is universal for given uni- 
versality class and boundary conditions. For liquid 4 He 
films it is believed that the superfluid order parameter 
vanishes at both surfaces of the film, implying Dirichlet 
boundary conditions Q . 

Unfortunately, a theoretical explanation of the strong 
dip and a determination of the scaling function is 
still lacking. In Ref. [l(|, this is stressed as the main theo- 
retical problem with respect to the explanation of the 4 He 
experiments. While field theoretical results 11, 12l [l3j 
are restricted to temperatures T >T C , Monte Carlo sim- 
ulations are only available for periodic boundary condi- 
tions until now [l4j, as only in this case the used stress 
tensor representation of the Casimir force is applicable. 
Recent attempts [ljl [l6[ to explain the strong dip within 
mean field theories only find qualitative agreement with 
the experiments, neglecting the non-critical contributions 
of Goldstone modes. 

In this letter, I present a direct calculation of the 
Casimir force using Monte Carlo simulations of the classi- 
cal XY model. This method requires the computation of 
the free energy of the system with high accuracy, which is 
a major challenge within Monte Carlo simulations. For- 
tunately, it turns out that the determination is greatly 
simplified by the fact that only the difference of two free 



energies is needed, which goes to zero exponentially fast 
above T c . The resulting finite-size scaling function is in 
excellent agreement with the experimental results l_ 
as well as with available theoretical predictions [^LOiri 
The Casimir force per unit area of a system with size 



Ljj 1 x L± (L|| — > oo ) is defined as [Til. IT 



pF Cas (T,L 1 _)=- 



0/cx(T, L A 



(1) 



where f c *{T, L±) denotes the excess free energy 

f ex (T,L x ) = f(T,L 1 _)-L ± f 00 (T) (2) 

of the system. Here f(T, L±) is the free energy per unit 
area of a film of thickness L±, measured in units of k^T, 
and foc{T) is the bulk free energy density. For large L± 
and near T c the Casimir force fulfills the scaling ansatz 

M 



/3F Ca ,(T,L ± )~Ll d tf(z) 



(3) 



with the universal finite-size scaling function ^(x) and 
the scaling variable 



l/u 



t>0 



L i 



l/u 



(4) 



Here I introduced the reduced temperature t = T/T c — 1 
and the bulk correlation length 



(i>0). 



(5) 



Note that the universal finite-size scaling function t?(x) 
depends on the boundary conditions in the Lj_-direction. 

We can directly calculate the Casimir force ([T]) by in- 
tegration of the internal energy as follows: Let us define 
the "internal Casimir force" 



pF ba (T,L x ) = - 



( du(T, Lj 



Uoo(T) 



with the internal energy per unit area in units of k^T 

,df(T,L±) 



(6) 



u(T,L±) 



-T- 



dT 



(7) 
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and the corresponding bulk density itoo (T) . The quantity 
du(T, L±)/dL± is directly accessible in Monte Carlo sim- 
ulations using u = (f3H)/L^, and the central difference 
quotient 



du{T,L ± ) u(T,L_l + 1)-u(T,Lx-1) 



dL ± 2 
The Casimir force is then obtained by integration 



(8) 



PF Cas (T,L ± ) = - f —l3F int (T,L ± ). (9) 
Jt t 

By Eqs. Q and ([9]), the internal Casimir force fulfills the 
scaling form 



(10) 



with the universal finite-size scaling function $'(x). Note 
that within the scaling regime, Eq. (fT0|) , the relative error 
of Eq. © is 0{L l / v ~ 2 ~ d ). 

We now consider the isotropic XY model on a simple 



cubic lattice of size L 



xLnxLt 



in three dimensions with 



periodic boundary conditions in the parallel directions. 
The Dirichlet boundary conditions in perpendicular di- 
rection are implemented by open boundary conditions, 
which are known to be equivalent at large length scales 
20l 21 1 , although alternative implementations are possi- 
ble 22 1 . The Hamiltonian reads 



H 



E 

(ij) 



Si * 



(11) 



where J > is the ferromagnetic exchange interaction, 
Si are 2-component unit vectors at site i, and the sum is 
restricted to nearest neighbors on the lattice. The simu- 
lations were performed for several system sizes with fixed 
aspect ratios p = L_\_/L\\ = 1:8 and 1:16 using the stan- 
dard Wolff cluster algorithm To calculate Eq. 
systems with thicknesses L'j_ = L± ± 1 at constant L\\ 
were simulated for every combination of Lii and Lj_. At 
least 10 5 Monte Carlo sweeps per data point were per- 
formed. 

To calculate the Casimir forces, Eqs. ([1]) and ©, it 
is necessary to have an expression for the bulk internal 
energy density itoo (T) . This is achieved using a combina- 
tion of direct simulations of a large cubic system (L = 96) 
with periodic boundary conditions and results of Cuc- 
chieri et al. 24jj . They determined the scaling behavior 



of the internal energy and specific heat of the XY model 
(fTTj) in the region \t\ < 0.015, where finite-size effects 
arise, using the scaling ansatz 



fc B Tuoo(T) = e ns + T c t 



— *" 
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1 



(12) 
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Figure 1: Results for the Casimir force using Eq. (114[) . for sys- 
tems with L± = 8 (green squares), L± = 12 (blue diamonds), 
and L± = 16 (red triangles), with aspect ratios p = 1:8 (open) 
and p = 1:16 (filled). The statistical error is of the order of 
the symbol size. The upper and lower insets are magnifi- 
cations around x = and around the minimum, respectively. 
Also shown are the Goldstone amplitude — £(3)/87T [2f (dashed 
line) and the field theoretical result pl| (cyan curve in upper 
inset), (color online) 



The critical indices of the considered XY model are fixed 
to the values v = 0.672(1), a = -0.017(3), u = 0.79(2), 
T c /J = 2.20183(1), and = 0.484(5) in the present let- 
ter, leading to the parameters e ns = —0.98841(3), C ns = 
22.03, A+ = 0.3790(8), A" = 0.3533(8), 4 = 0.015(1), 
c{ = 0.109(2), c+ = -0.041(3), and = 0.211(4) Q. 
For \t\ > 0.015 finite-size effects are negligible; note that 
at t = 0.015 the correlation length, Eq. ([5]), has the value 
(0.015) « 8.1, which is sufficiently small with respect 
to L = 96. While for periodic cubic systems the scaling 
corrections are moderate, systems with broken transla- 
tional invariance and aspect ratios p< 1 show strong cor- 
rections to scaling. An analysis of usual thermodynamic 
quantities like the magnetic susceptibility x(T, L±) and 
the Binder cumulant U (T, L±) shows that it is necessary 
to use a modified scaling variable x (Eq. (0|) with Wegner 
corrections [25( of the form 



l/u 

^ = Mtt) (i + g u Li u ), 



£ + 



(13) 



while the y-direction has rather small corrections for sys- 
tems with constant p. However, for the numerical deriva- 
tive with respect to L± in (J6]) it is necessary to combine 
data of systems with L' ± = L± ± 1, leading to systems 
with different aspect ratio p' ^ p. This and the expected 
uncertainty of the numerical derivative itself introduces 
a scaling correction of the order (1 + giLJ 1 ) in the y- 
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Figure 2: Universal finite-size scaling function #(a;) of the Casimir force for systems with L± = 8 (green squares), L±_ = 12 (blue 
diamonds), and L± — 16 (red triangles), with aspect ratios p = 1:8 (open) and p = 1:16 (filled). The insets are magnifications 
of the respective regions. The results are compared (left) to the experimental data of Garcia and Chan [y, Cap. 1] (•), as well 
as with the results (right) of Ganshin et al. (A : 340A, ♦ : 285A, T : 238A; solid lines are guides to the eye). Also shown 
are the Goldstone amplitude — £(3)/8n || (dashed line), the value — llC(3)/327r including surface fluctuations proposed in 0] 
(dotted line), and the field theoretical result [l^] (cyan curve in left upper inset), (color online) 



direction, leading to the final scaling ansatz 

l3F Cas (T,L ± ) ^ L- d (l + g.L- 1 )- 1 ^) (14) 

with x from Eq. fl3]). 

The results for the Casimir force are shown in Figure [1] 
for six system sizes with L± G {8, 12, 16}, each with as- 
pect ratio p = 1:8 and p = 1:16. It should be emphasized 
that the only fit parameters are the corrections to scaling 
amplitudes, g^ = 2.0(1) and g\ = 5.5(2), which are ad- 
justed until the numerical data collapse onto a single 
curve, and that there is no free factor in neither x- nor 
y-direction. We can identify a slight dependence on the 
aspect ratio p in both directions. As the corrections due 
to the finite p are known to scale approximately with p 2 
[27l ] , a full data collapse can be achieved by adding a fac- 
tor (1 + r\p 2 ) to the x-axis and a factor (1 + r 2 p 2 ) to the 
y-axis, with r\ = 4(1) and r 2 = 10(1). Note that these 
corrections mainly shift the curves for p = 1:8, while the 
p = 1:16 curves are virtually unchanged within the er- 
ror bars. The resulting scaling function is depicted in 
Figure O together with the results of Garcia and Chan 
(left) as well as Ganshin et al. [8j (right). In the 
first case, only data for the thickest film with d = 423 A 
are shown, which are regarded to have the highest qual- 
ity [28}, while the thinner films showed deviations in y- 
direction, presumably due to surface roughness @, @] . In 
order to compare the experimental data quantitatively 
with the present results, they are made dimensionless in 
x-direction using the measured correlation length ampli- 
tude ^ = 1.432 A of 4 He at T A [29], leading to a factor 



(^y 1 ^ = o^eA- 1 /^. 

We find an excellent agreement within the error bars 
with both measurements for x > —8. The universal am- 
plitude of -d(x) at the minimum is ^(a; m i n ) = —1.35(3) 
at x mnl = —5.3(1). These values agree with the values 
a^min = —5.4(1) of Garcia and Chan 0] and i?(a; m in) = 
— 1.30(3) at x m in = —5.7(5) of Ganshin et al. Q. It turns 
out that the overall agreement with is even better than 
with which might be attributed to the smaller fluc- 
tuations in y-direction in the data of 0], mainly visible 
below x w —10, and to the five times smaller error es- 
timate in x-direction, clearly visible in the insets and at 
x = 0. For x < — 8 we see an enhancement of the mea- 
sured Casimir force not present in the calculated scaling 
function. This onset is weaker in the left figure. A possi- 
ble explanation is the occurrence of surface fluctuations 
below this temperature, as proposed in Q. At the critical 
point (x = 0) we find i?(0) = -0.047(2) [-0.059(2)] for 
p = 1:8 [1:16], which gives 0(0) = -0.062(5) for p -> 0. 
The resulting Casimir amplitude A = —0.031(3) agrees 
well with the estimate A = -0.03 from [!,[27j]. For x > 1 
the results nicely lie on the scaling function calculated 
by Krech and Dietrich [12 ] using renormalization group 
theory. Here r\ — 30 was used, which again mainly shifts 
the data for p — 1:8, see upper inset in Figure [2jlcft) . 
The quality of the used method, especially of the numer- 
ical integration, is demonstrated by the convergence of 
the calculated scaling function to the low temperature 
Goldstone value i?Goldstone = — C(3)/87r [2] for x — > —00 
(dashed line in Figures [1] and \2§ ■ 
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Figure 3: Universal finite-size scaling function $'(x) of the in- 
ternal Casimir force, Eq. ([6]), determined from eight different 
system sizes with aspect ratios p = 1:8 and 1:16. The exper- 
imental results are obtained by binning (with Ax = 0.5) and 
numerical differentiation of the data from Garcia and Chan 
0, Cap. 1]. (color online) 



Figure [3] shows the results for the scaling function 
i)'(x) of the internal Casimir force, Eq. ((6]). The scal- 
ing plot contains data from eight system sizes with L±_ G 
{4, 8, 12, 16}, each with aspect ratio p = 1:8 and p — 1:16. 
While the aspect ratio correction in y-direction becomes 
f*3 = V2 — Ti, the scaling correction gi cannot be ex- 
pressed through and <?i, as the corresponding L±- 
exponents are different. The effective calculated value 
at L± w 10 is 92 = 1-7, which is modified to 32 = 2.0(2) 
to get the best data collapse. The results are compared 
to a numerical differentiation of the experimental data 
of Garcia and Chan 0, Cap. 1]. The results of Ganshin 
et al. 8] are not shown due to large fluctuations of the 
numerical derivative. The data collapse and the agree- 
ment with the experimental data is very convincing, also 
showing the small influence of statistical errors in the 
simulations. Only in the interval —9 < x < —5 we see 
higher order corrections to scaling, which are believed 
to stem from uncertainties in the numerical derivative, 
Eq. (jHJ. However, these corrections only have small in- 
fluence on the integrated Casimir force. Further work is 
necessary to clarify this behavior. Note that at the min- 
imum of $(x) we have $'(x m i n ) = 0, which implies that 
\im L± _ 00 F int (T min (L ± ),L ± ) = and 



du(T, Lj 



dL± 



Moo C^min (£_!_)), 



(15) 



T=T min (Zx) 



i.e. at this temperature (and large L±) the change in 
internal energy with L± equals the corresponding bulk 
value. The connection (1151) between the minimum of the 



Casimir force and the internal energy shows that the shift 
of the minimum to negative a; is a direct consequence 
of the strong shift in T C (L±) in systems with Dirichlet 
boundary conditions. This shift is also present in the 
exactly solvable two dimensional Ising model [l8j |. 

In summary, I determined the universal finite-size scal- 
ing function i3(x) of the Casimir force within the XY 
universality class with Dirichlet boundary conditions us- 
ing Monte Carlo simulations. For sufficiently small as- 
pect ratio p = 1:16 the results are in excellent agree- 
ment with the experimental results on 4 He by Garcia 
and Chan Q, and by Ganshin et al. [&], as well as with 
theoretical calculations for T > T c by Krech and Dietrich 
ll|,ll2j. The universal function has a deep minimum 
at x min = —5.3(1), with $(:c m i n ) = —1.35(3). The results 
are in conformity with the assumption that the order pa- 
rameter in 4 He asymptotically obeys Dirichlet boundary 
conditions. The method proposed in this letter has also 
been applied to systems with periodic boundary condi- 
tions [30(| , where a similar good agreement with avail- 
able results 13, 14 1 is obtained. The application to other 
boundary conditions as well as to Ising and Heisenberg 
models is straightforward. 

Note added in proof. Recently, Vasilyev et al. pre- 
sented an alternative method to calculate $(x) using 
Monte Carlo simulations [3l|. 
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